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Quantum properties of optical modes are typically assessed by observing their photon statistics or the distribu¬ 
tion of their quadratures. Both particle- and wave-like behaviours deliver important information, and each may 
be used as a resource in quantum-enhanced technologies. Weak-field homodyne detection provides a scheme 
which combines the wave- and particle-like descriptions. Here we show that it is possible to observe a wave-like 
property such as the optical coherence across Fock basis states in the detection statistics derived from discrete 
photon counting. We experimentally demonstrate these correlations using two weak-field homodyne detectors 
on each mode of two classes of two-mode entangled states. Furthermore, we theoretically describe the response 
of weak-field homodyne detection on a two-mode squeezed state in the context of generalised Bell inequalities. 
Our work demonstrates the potential of this technique as a tool for hybrid continuous/discrete-variable protocols 
on a phenomenon that explicitly combines both approaches. 


Quantum correlations play a central role in our understanding of fundamental quantum physics and represent a key re¬ 
source for quantum technologies [1-4]. Progress in quantum information science has followed an increasingly thorough 
understanding of how such correlations manifest themselves, and how they can be successfully generated, manipulated and 
characterised [5]. In quantum optical systems, these correlations appear in either properties of the fields, such as quadrature 
entanglement, or at the level of individual quanta, i.e., photons. Access and control over such correlations are key to appli¬ 
cations in quantum information science; however, how the complementary mode- and particle-correlations precisely act as a 
resource is still a subject of debate [6]. The capability of studying the correlations in a regime traversing mode and particle 
aspects is thus necessary for clarifying the origin of quantum enhancement. 

When a coherent superposition of many photons occupies a single mode, a wave-like description of the quantum state in 
terms of continuous variables (i.e., the values of the quadratures) of the electromagnetic field is the standard approach [7]. The 
canonical technique for measuring such light fields is strong-field homodyne detection, which directly probes the quadratures 
of a field and can provide a full reconstruction of its quantum state [8]. On the other hand, particle-like properties can be 
directly accessed with a range of photodetectors, a notable example being photon-number-resolving detectors (PNRDs) [9- 
12], and suitable techniques for the reconstruction of the photon statistics [13, 14]. However, such photon counters are 
intrinsically insensitive to phase, and thus cannot access any coherence between modes. Weak-field homodyne (WFH) 
has been proposed as a measurement technique bridging wave-like and particle-like descriptions [15-21]. As in standard 
homodyne, a local oscillator (LO) provides phase sensitivity, while the photon statistics are accessed by number-resolving 
detectors. The main difference is that WFH makes use of a classical reference whose mean photon number is of a similar 
order as that of the probed signal; this allows for the combination of the homodyne technique with PNRDs based on photon 
counting modules. 

In this paper, we employ WFH detection to investigate coherence between different photon-number basis states (Fock 
layers) across two-mode entangled states. Our detection scheme accesses this manifestation of optical coherence directly, 
without the need for resource-intensive full state tomography. We demonstrate the oscillations of an array of multi-photon 
coincidence counts when a split single-photon state (SSPS) and a two-mode squeezed state (TMSS) are interfered with a 
weak local oscillator. Our experiment can be regarded as a first step towards the quantitative study of the nonlocal properties 
of multi-mode quantum states with multi-outcome measurements by non-Gaussian detectors. To this end, we theoretically 
study a violation of Bell inequalities that can be achieved for low numbers of detected photons under ideal conditions. 


RESULTS 

Experiment. The weak-field homodyne setup is depicted in Fig. la. A signal is mixed on a beam splitter with a local 
oscillator, described by the coherent state |o:e^), which establishes a phase reference. One or both output modes from 
the beam splitter are then detected using photon-number-resolving detectors. The method is termed “weak-field” since the 
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FIG. 1 : Layout of the weak-fieid homodyne (WFH) detector, (a) The general scheme relies on the interference between a signal and a 
local oscillator (LO) of similar intensity for different settings of the phase 0. Either one or both outputs are measured with photon-number- 
resolving detectors (PNRDs). (b) Our experimental implementation adopts a collinear configuration in which the signal |^) and the LO 
have orthogonal polarisations. For this reason, the beam splitter (BS) in (a) is replaced by a half-wave plate (HWP) and a polarising beam 
splitter (PBS) which realise interference in polarisation. The phase setting of the LO is controlled by means of a geometric phase rotator 
(GPR), which consists of a quarter-wave plate (QWP), a half-wave plate and a second QWP. The rotation of the HWP (with both QWP 
fixed at 45°) applies a phase shift to the coherent state. A time-multiplexed detector (TMD) records clicks on the transmitted output mode. 


intensity of the LO is comparable to that of the signal. This differs from strong-field homodyne in which the coherent state is 
many orders of magnitude more intense than the signal, and the outputs are detected by linear photodiodes. 

To understand how WFH can detect coherence across Fock layers, consider the detection of m photons. These can originate 
from any Fock term \k) in the signal together with m—k photons from the LO; if the complex coefficients of the terms {| k)} 
have well-defined relative phases, then these result in a modulation of the detection probability p(m) as the phase of the local 
oscillator is varied. For a two-mode state, the relevant phases are the ones between joint detection events coming from terms 
in the signal of the form |&i)|& 2 )- Consequently, coherence across such events can be observed in the modulation of joint 
detection probabilities. 

We use the layout shown in Fig. lb to study the coherence between Fock layers across the two modes of a split single 
photon and a squeezed vacuum state. To maximise the passive phase stability of our setup, we adopt a compact design in 
which the LO and the signal are collinear and occupy two orthogonal polarisations of a single spatial mode. Our PNRDs are 
time-multiplexed detectors (TMDs) that split an incoming beam into two spatial and four temporal modes, thereby resolving 
up to eight photons probabilistically using two avalanche photodiodes [10, 22]. Specifically, the TMDs allow us to decompose 
the intensity of the interference patterns resulting from each output mode into its constituent photon components. In this way 
we can probe pair-wise correlations between individual Fock layers as we build up a joint detection statistics matrix, every 
row and column representing the number of clicks in each detected mode. The click statistics gives us access to higher-order 
Fock states (k > 1), although this detection scheme is not fully equivalent to a number-resolving detector [23]. 




FIG. 2: Probing correlations across optical modes with WFH. (a) A single photon is produced by a heralded source based on parametric 
downconversion (PDC). We then split the single-photon state in order to generate |T^)ssps across two separate spatial modes (A and B), 
both probed by the WFH detectors described in Fig. lb. (b) The same downconversion source produces the squeezed vacuum state, whose 
two modes are also probed by weak-field homodyne detectors. 
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A split single-photon state may be written as 

|^)ssps = —^ (|0) a |1)b + |1)a|0)b), (1) 

which describes a single photon in a coherent superposition of modes A and B (see Fig. 2a and Supplementary Fig. 1 for 
further details on the experimental layout). The results of our investigation for the SSPS are reported in Fig. 3. 

The red dots show the experimental probabilities P(ra, m') of the joint detection of m clicks in mode A and m' clicks 
in mode B. Each P(ra,ra') term is a function of the difference between the phase settings of the weak local oscillators, 
since a single photon has no phase per se [24, 25]. Our count rates are such that we only consider events 
where m,m' < 2; detection outcomes greater than this level are negligible. The blue lines are predictions from a model 
which accounts for the following imperfections in our experiment: non-unit efficiency of the detectors; modulation of the 
reflectivities of the beam splitters preceding the time-multiplexed detectors when varying the LO phases (due to the geometric 
phase rotators depicted in Fig. lb - see Methods and Supplementary Fig. 2). We also include imperfections in the production 
of the single-photon states [26], and we thus write 

Po = wo|0)<0| + iui|l)(l| + (1 - u; 0 - wi)|2)<2|, (2) 

where wo and w\ are coefficients taken between 0 and 1 which weight the zero- and one-photon contributions to the input 
state. We experimentally determine these parameters from the photon statistics of the initial state (see Supplementary Note 1 
and Supplementary Note 2). 
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FIG. 3: Joint counting statistics for the SSPS. Correlations between the responses of two WFH detectors as a function of the phase 
difference = 0^/4 (where is determined by the settings of the half-wave plates in the phase rotators). We recall that 

P(m, m') is associated to the joint detection probability of m photons on mode A and m' photons on mode B. The red dots are measured 
probabilities, with uncertainties determined by a Monte Carlo routine under the assumption of a Poissonian distribution for the raw counts. 
The blue theoretical curves are obtained from a model which includes the main imperfections in the setup. The dashed curves correspond 
to an input state with the experimentally determined weights u;o=0.161 and w\ =0.669. The uncertainties cro=0.011 on wo and g\ =0.029 
on w i are estimated by a Monte Carlo routine; the dotted curves show the results given by the theoretical model with (wo — cjo, w± — cr i) 
and (wo+ao, wi + cri), respectively. The quantum efficiencies of the detectors, tja — 0.072 and tjb — 0.064, are experimentally estimated 
with the Klyshko method [27, 28]. The intensities of the two local oscillators are |qa|= 0.510 and |qb|= 0.585. 

The experimental curves show oscillations in the coincidence counts: these are evident in the P(1,1),P(2,1) and P(l, 2) 
terms, as also predicted by the theory. What is most striking about the observed oscillatory behaviour is the fact that it is 
displayed by terms P(m, m') for which m + m! > 1, yet it is determined by the coherence between the terms |0)a|1)b an d 
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| 1)a|0)b (see Eq.(l)) which do not contain more than one photon each. Ideally, any additional photon detected must therefore 
come from the local oscillators, hence the SSPS is responsible for the coherent oscillations observed in the considered joint 
detection probabilities. In practice, we observe good qualitative agreement between the experimental data and our theoretical 
description. Indeed, we are able to account for all the main features of a 3-by-3 array of multi-photon counts with one model 
that has no free parameters: detection efficiencies and the weights in Eq.(2) are experimentally determined. We attribute the 
residual discrepancies to two main factors: imperfect mode matching between the LOs and the signal modes, and variations 
of the laser power during the few-hour long acquisition. 

An analogous experiment is conducted on a two-mode squeezed state (see Fig. 2b and Supplementary Fig. 2 for further 
details on the experimental layout). This state is archetypal of a quantum resource with well-known correlations between its 
Fock layers. The expression for a TMSS reads 


+oo 

|^)tmss = yr^£A>) A |n>B, (3) 

n =0 

where the real squeezing parameter | A| governs the photon distribution across the photon-number basis states. Here we pump 
our source with moderate parametric gain in order to generate significant higher-order Fock layers in the two-mode state [29] . 
The expected phase dependence of {P(ra, m')} for a two-mode squeezed state was shown to be 0(+)=0 a +0b [30]. The 
TMSS has a phase dependence arising from that of the pump; in a conventional strong homodyne setup this would define 
which quadratures are squeezed. In our case, the same effect is manifested in the phase dependence of the click patterns; 
hence, appropriate phase locking was necessary (see Methods). 
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FIG. 4: Joint counting statistics for the TMSS. Correlations between the responses of two WFH detectors as a function of the phase 
sum = (j >^/4 (where is determined by the settings of the half-wave plates in the phase rotators). We recall that P(m, m) is 
associated to the joint detection probability of m photons on mode A and m photons on mode B. The red dots are measured probabilities, 
with uncertainties determined by a Monte Carlo routine under the assumption of a Poissonian distribution for the raw counts. The blue 
theoretical curves are obtained from a model which includes the main imperfections in the setup. The dashed curves correspond to an input 
state with the experimentally determined squeezing parameter | A|=0.295 and weight p=0.04 for the additional noise term. The uncertainty 
cr=0.016 on | A| is estimated by a Monte Carlo routine; the dotted curves show the results given by the theoretical model with | A| =b a. The 
experimental quantum efficiencies of the detectors are t]a = 0.132 and t]b = 0.155 (estimated with the Klyshko method). The intensities 
of the two local oscillators are |a A | = 0.365 and |qb|= 0.347. 


The results obtained for the two-mode squeezed state are illustrated in Fig. 4. Once again, we observe good qualitative 
agreement between our data and the predictions from a theoretical model that includes the same imperfections as for the 
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single-photon case. Here the input state can be modelled as 

p = (1 - p)|^)tmss(^| +p|0,1)(0,1|, (4) 

where the extra term is a first-order approximation of noise in a squeezed thermal state. This asymmetry across the modes 
is justified by the experimentally recorded g^ for the parametric downconversion source: we find g^ = 1.98 ± 0.04 and 
g = 1.92 dz 0.04 for the two marginals, to be compared with g^ = 2 for an ideal thermal state. The lower value of 
g motivates the addition of the noise term on output mode B. The quoted values of g^ also suggest that additional (and 
undesired) Schmidt modes might be responsible for the presence of photons in modes correlated to undetected modes [31]. 
As for the case of the SSPS, imperfect mode matching between the LOs and the signal modes and variations of the laser 
power during the data acquisition are recognised as the main causes of the residual departures of the experimental curves 
from the expected behaviour in Fig. 4. 

As a general remark applying to both classes of studied states, we note that our numerical models depend strongly on the 
detection efficiencies (see Supplementary Fig. 4 and Supplementary Fig. 5). In this sense, the more pronounced discrepancies 
in Fig. 3 and Fig. 4, such as that in the P(0, 2) term, may be due to noise affecting our estimation of these parameters. More 
details on the theoretical model for both the SSPS and the TMSS can be found in Supplementary Note 1 and Supplementary 
Note 2. 

Generalised Bell inequalities for WFH detection. Correlations such as those revealed in our experiment need to be 
understood in relation to canonical criteria for non-classicality, for instance the violation of a Bell inequality, to assess their 
role as possible quantum resources [32]. The ability to access and discriminate higher photon numbers leads us to refer to 
generalised, higher-dimensional Bell inequalities [33], which we study theoretically in the context of our experiment with a 
TMSS. 

To this end, let us consider the scenario depicted in Fig. 5. Each mode of a TMSS is analysed by means of a WFH 
detector; here all four outputs are monitored by perfect PNRDs. We take into account Fock layers which lead to the detection 
of M photons on each side: the detected photons are split among the two output ports of the PBSs (i.e., the photons are 
either transmitted or reflected) according to the convention indicated in Fig. 5a. Specifically, we consider measurements with 
D=M+ 1 possible results, distinguished according to the number of photons T detected on the transmitted arm of each PBS. 
We are thus interested in the probabilities associated with joint detection events comprising outcomes (on each output arm, 
on both sides) differing by e (the outcomes being taken modulo D). {a, a'} and {/?,/?'} denote the LO settings on side A 
and side B , respectively. The relevant probabilities are combined into the expression 

m-V 2e\ r 1 

Im = E \l-^j[P(^/3,e) + P(a , ,/3,-e-l) + P(a , ,l3 , ,e) + P(a,/3 , ,-e)\ ^ 

— ^P(a, /?, —e — 1) + P(a', /?, e) + P(ol , /3', —e — 1) + P(oi, /3' , 6 + 1) , 

where the local realistic bound is \Im\ < 2 [33]. We note that the case M= 1 corresponds to the standard CHSH inequality 
which was experimentally tested in [17]. Additional details on how to compute Im for the specific layout that we consider 
are provided in Supplementary Note 3, Supplementary Fig. 6 and Supplementary Fig. 7. 

We run a numerical search for values of M ranging from 2 to 8 in order to find the set of parameters {A, a, a/, f3, /?'} 
which determine the highest violation of these generalised Bell inequalities. Our results are shown in Fig. 5: no violation of 
Eq.(5) is found beyond M= 5. This behaviour comes from the fact that the maximal violation of Eq.(5) (i.e., the attainment 
of the maximum value allowed by quantum mechanics) relies on a particular structure for the entangled state [34]. Notably, 
a two-mode squeezed state is specified solely by the squeezing parameter A: this means that there are not enough degrees 
of freedom to tune in order to obtain the required form for the input state. This restriction becomes more severe as the 
dimensionality of the system increases, to the point that no violation of local realism can be inferred despite the use of an 
entangled resource such as a TMSS. Hence the decrease in the violation stems from the Gaussian character of the two-mode 
squeezed state, which fixes the functional shape of the oscillation curves simultaneously for all Fock layers. 


DISCUSSION 


We have shown how WFH can be used to deconstruct phase-dependent measurements on two-mode entangled states into 
their constituent Fock layers. The ability to operate devices at the interface of wave and particle regimes opens up new 
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FIG. 5: Generalised Bell inequalities with WFH. (a) In the ideal case, the two modes of a TMSS are interfered with two weak local 
oscillators; PNRDs monitor all four outputs. We fix the number of detected photons to M on each side, and label the various detection 
events recorded by the PNRDs according to how the M particles distribute themselves on the two output modes: r, M — T, T + e and 
M—T— e. (b) The plot shows the values taken by the quantity Im (see Eq.(5)) when the number of detected photons M on each side varies 
between 2 and 8. The bound for Im predicted by local realism is 2 for all values of M, hence we see that no violation of a generalised Bell 
inequality is possible beyond M — 5. 


possibilities for quantum information processing [35]. Thus, our work provides a new insight on such resources within the 
broader investigations on hybrid continuous/discrete-variable coding [36^0]. 

As a complement, we have also studied the theory of the violation of generalised Bell inequalities using an ideal WFH 
setup. These tests shed light on how the transition from non-Gaussian to Gaussian measurements occurs. Such transition is 
achieved in the context of WFH detection as an increase in the LO strength. In fact, this study may also be interpreted in the 
more general framework of Gaussian vs non-Gaussian measurements, where it is well-known that the outcomes of a Gaussian 
measurement on a Gaussian quantum state can be explained by a local realistic model. For this reason, strong-field homodyne 
detection on a TMSS cannot be used to violate a Bell inequality [41, 42]. On the other hand, WFH is an example of non- 
Gaussian measurement, as attested by its Wigner function. Consequently, we expect that the outcomes of WFH detection on 
an entangled Gaussian state cannot be described by a local realistic model [17, 43, 44]; however, this breaks for moderately 
high photon numbers. 

On the experimental side, further developments for the observation of the violation of higher-dimensional Bell inequalities 
demand detectors with higher quantum efficiency [45]. These are not only necessary for achieving significant counting 
statistics, but also to prevent one Fock layer to be affected from higher-order contributions. Encouraging results have been 
obtained in this direction with cryogenic detectors. 


METHODS 

Source of quantum states. A pulsed Ti:Sapphire laser (repetition rate 256kHz, central wavelength AxiSa = 830nm and 
bandwidth A AxiSa — 30nm) is doubled in a nonlinear crystal (/5-barium borate, BBO). The second-harmonic beam (Auv = 
415nm) pumps a type-II parametric downconversion process in a nonlinear crystal (potassium dihydrogen phosphate, KDP) 
in order to produce a two-mode vacuum-squeezed state. The source is designed to generate spectrally uncorrelated modes, 
based on group velocity mismatch inside the birefringent medium [29] . Daughter photons have orthogonal polarisations and 
different spectral widths (AAy ~ 12nm, AAr — 6nm). The same source generates the split single-photon state and the 
two-mode squeezed state, depending on the ultra-violet pump beam brightness. 

Weak-field homodyne detection. The detection system that we adopt in our experiment is realised with a collinear 
geometry in order to ensure passive phase stability. Each mode in Eq. (1) and (3) is superposed with a local oscillator of 
orthogonal polarisation at a polarising beam splitter (PBS). This delivers a common spatial mode at the output but orthogonal 
polarisations for the signal to be probed and the weak coherent beam. Interference is then realised by a half-wave plate and an 
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additional PBS. The relative phase between the horizontal and vertical polarisations can be adjusted by means of a geometric 
phase rotator (GPR - see Fig. 1), which is composed of a quarter-wave plate (QWP), a half-wave plate and a second QWP. 
The axis of the first QWP is set to 45° in order to transform the input linear polarisations into circular ones. The rotation of 
the HWP by an angle 0/4 results in a phase shift equal to 0 between left- and right-circular polarisations. The initial linear 
polarisations are then recovered by setting the second QWP to 45°. The successful calibration of the full device relies on the 
characterisation of each element (including the PBS), particularly of the QWPs. When the GPR is correctly calibrated power 
fluctuations around 3% are recorded, while imperfections in the calibration of one of its constituents can lead to fluctuations 
in power above 10% (see Supplementary Fig. 2). Finally, the output state is analysed with a time-multiplexed detector: this 
consists of two fibre-based cascaded Mach-Zehnder interferometers that split the incoming state over two spatial modes and 
four distinct time bins. Time-resolved clicks from avalanche photodiodes monitoring the two transmitted modes are thus 
registered. 

Active phase stabilisation. To actively lock the phase set by the GPRs, we use an ancillary laser beam (from a continuous- 
wave HeNe laser, A' = 633nm) which back-propagates through the interferometer. The classical interference pattern that 
we obtain when the ancillary beam reproduces correctly the signal and LO optical paths constitutes the signal recorded by a 
photodiode connected to a PID device (SRS SIM960 Analog PID Controller). 
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SUPPLEMENTARY INFORMATION 
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Supplementary Figure 1: Experimental layout for the study of optical coherence via weak-field homodyne across the modes of a 
split single-photon state (SSPS). A blue beam (Auv = 415nm) interacts with a Type-II potassium dihydrogen phosphate birefringent 
crystal to produce frequency-degenerate downconverted photons with orthogonal polarisations at 830nm. The ultraviolet beam is obtained 
via second-harmonic generation on a Type-I /3-barium borate crystal pumped by the radiation coming from a pulsed femtosecond 
Ti:Sapphire laser system (Coherent Mira Seed - RegA, central wavelength AxiSa = 830nm, repetition rate 256kHz). The vertically 
polarised downconversion mode is detected by an avalanche photodiode and thus heralds the presence of a signal on the twin mode, which 
is mixed with a local oscillator on a polarising beam splitter (PBS). The signal and coherent state, now collinear, are then interfered in 
polarisation on two PBSs preceded by geometric phase rotators which allow us to vary the phase of the local oscillators independently on 
Alice’s and Bob’s side. The two transmitted output modes are monitored by two time-multiplexed detectors, while the reflected ones are 
discarded. A field-programmable gate array processes the recorded coincidence and single counts. BBO: /3-barium borate crystal - DM: 
dichroic mirror - WP: wave plate (yellow, quarter-wave plate; red, half-wave plate) - PBS: polarising beam splitter - KDP: potassium 
dihydrogen phosphate crystal - APD: avalanche photodiode - GPR: geometric phase rotator - TMD: time-multiplexed detector. 
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Supplementary Figure 2: Effect of an imperfect geometric phase rotator (GPR) on the marginal distributions for detection events 
up to two photons, (a) Theoretical marginal distributions for the split single-photon state, (b) Experimental marginal distributions for the 
SSPS when the GPR on Bob’s side is moved, (c) Typical calibration curve for a geometric phase rotator. Theoretically, no term in the 
marginal counting statistics should oscillate when varying the phase of the local oscillator (see plot (a)). However, the experimental data 
(see graph (b)) reveals the presence of oscillations in the recorded counts on the mode where the geometric phase rotator is active. 
Analogous results are found for the TMSS. These plots suggest that a non-ideal GPR can cause spurious oscillations in the marginals, that 
is, the distributions of detection events when we trace over one mode of the full two-mode output state. The calibration curve (see plot (c)) 
is obtained by interfering the vacuum with a local oscillator and recording counts on one detector as the GPR is moved. We note that the 
visibility of the oscillations in the LO counts is about 10% for the 1-click term. If we compare this value with the visibilities of the 
oscillations in plot (b), we see that it falls between the 1-click term (V ~ 6%) and the 2-click term (V > 10% - although in this case the 
absolute values are of the order of 10 -4 counts • s -1 ) associated to the experimental marginal distributions. These observations support 
the addition of a modulation of the reflectivity in our theoretical models for the split single photon and the two-mode squeezed state. 
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Supplementary Figure 3: Experimental layout for the study of optical coherence via weak-field homodyne across the modes of a 
two-mode squeezed state (TMSS). The setup for the generation of two-mode squeezed vacuum states exhibits a few modifications with 
respect to that illustrated in Supplementary Fig 1. In this case we interfere both downconversion modes with a local oscillator (LO); 
indistinguishability in polarisation is achieved thanks to a half-wave plate and a polarising beam splitter on each output mode. Once again, 
the geometric phase rotators allow us to vary the phases of the local oscillators. This layout requires active phase stabilisation to ensure 
that the phase difference between each initial LO phase and the phase of the two-mode squeezed state is kept fixed. Our active phase 
locking relies on an ancillary laser beam from a continuous-wave HeNe laser (A ; = 633nm) which back-propagates through the 
interferometer. The signal produced by the classical interference pattern - which is obtained when the ancillary beam reproduces correctly 
both optical paths in the interferometer - is recorded by a photodiode connected to a PID device (SRS SIM960 Analog PID Controller) 
whose electronic signal is then fed to a piezo delay stage in the LO path (indicated by the gray dashed box in the diagram). BBO: 
/3-barium borate crystal - PBS: polarising beam splitter - HWP: half-wave plate - KDP: potassium dihydrogen phosphate crystal - GPR: 
geometric phase rotator - HeNe: helium-neon diode laser - PD: photodiode - PID: Proportional-Integral-Derivative device. 
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Supplementary Figure 4: Theoretical joint counting statistics for different detection efficiencies in the case of a split single-photon 
state. To illustrate the effect of imperfect detection on the observed nonclassical oscillations, we plot the array of multi-photon 
coincidence counts while varying the values of the detection efficiencies. Green curves correspond to 77 = 0.05, cyan to 77 = 0.1 and blue 
to 77 s= 0 . 2 ; we assume the same efficiencies on both time-multiplexed detectors and use the experimentally determined values of the 
coefficients wo and w\. We see how the parameter 77 changes both the level of detected signal and, most importantly, the visibility of the 
oscillations. 
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Supplementary Figure 5: Theoretical joint counting statistics for different detection efficiencies in the case of a two-mode 
squeezed state. We illustrate the effect of imperfect detection on the observed oscillations by varying the values of the detection 
efficiencies. Green curves correspond to 77 = 0.05, cyan to 77 = 0.1 and blue to 77 = 0.2; we assume the same efficiencies on both TMDs 
and use the experimentally determined values of the squeezing parameter | A| and weight p. As in the case of the SSPS, the parameter 77 
determines dramatic changes in the visibility of the nonclassical oscillations. 
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Supplementary Figure 6: Study of Gaussian entanglement with weak-field homodyne (WFH) detection. This is the ideal layout for 
the investigation of the nonlocal character of Gaussian entanglement via WFH. A two-mode squeezed state, an archetypal Gaussian 
resource, is interfered with a weak coherent state on two PBSs; photon-number-resolving detectors Di, Ei, D 2 , E 2 monitor all four 
outputs (i.e., we do not discard the reflected output modes here). If we fix the number M of photons detected by each WFH detector, we 
can label the recorded detection events according to how M particles distribute themselves on the transmitted and reflected output modes: 
namely, these are the quantities T, M — T , T +e and M — T — e shown in the diagram. 
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Supplementary Figure 7: Probabilities of four-photon detection events as functions of the phase of the local oscillator. The graphs 
show the behaviour of the expectation values of the projection operators on the considered four-photon states, i.e., the set {( P^k ')} for 
k, k! < 2 (see Supplementary Note 3): here we take into account the multiplicative factor yj 1 — |A| 2 e - ^ neglected in the 

expression for |T0 O ut (see Eq.(23)). The phase on the horizontal axis is that of either one of the local oscillators, 0 a or 0&; we set 0 = 0 for 
the two-mode squeezed state and 0& = 0 a + tt/ 4, relabelling <j) a — 0. Physically plausible values - under the assumption of ideal 
detectors, i.e., rj = 1 for all PNRDs - are assigned to the other variables so that the considered functions depend solely on one phase 
parameter: we have |A| = 0.3 and |a| = \/3\ = 0.131. We also assume the PBSs to be balanced, therefore U — rj = 1/V% for i, j = 1,2. 
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Supplementary Note 1: Detection probabilities for the weak-field homodyne 


We describe the expected response of a weak-field homodyne (WFH) detector to a generic quantum state. For this purpose, 
we study the interference of a Fock state | n) with a coherent state | a) on a beam splitter (BS) with transmittivity t and 
reflectivity r = \/l — t 2 . We follow the asymmetric convention for the beam splitter input-output relations, so that the modes 
are transformed as aj n faj ut + rb\ nV and bj n tbj ut — raj ut . From this, it follows that the operators which generate 
coherent and Fock states transform as D in (a) D oat2 {ta)D outl (-ra ) and (aj,) n -s- YJk =o (l)t n ~ k r k (al utl ) n ~ k {al at2 ) k , 
respectively. Therefore, given the initial state 

|^) in = |a)|n>=£> in (a)(aj 1 ) n |0)|0) 1 (6) 

this undergoes the above-described beam splitter transformations and becomes 

Wout = £ Wf"- fc r fe J D 0Ut2 (to) J D 0Utl (-m)(dJ utl ) n - fe (aL 2 ) fc |0)|0>. (7) 

k=0 ' ' 

We introduce the short-hand notation |a; n) := D(a) (al) n |0) |0) for the displaced Fock state; we note that the norm of \a; n) 
is n\. For a given displacement a , these states form an orthogonal set, (a; n\a\ m) = \fn\ S nm , as the displacement preserves 
the scalar product. This notation allows us to express the total output state in Eq.(7) as 

Wout = V ( T f\ t n ~ k r k \ta; n — k)\-ra;k). (8) 

k =o ^ ' 

As a consequence, an input state whose expression in the Fock basis is c n\ n ) transforms as 

+oo + 00 / \ 

|^)out = ^2^2 c n ( r ^jt n ~ k r k \ta;n - k)\ -ra;k), (9) 

k =0 n=k ' 

where we have inverted the order of the summations over n and k. Our experimental detection scheme does not record the 
outcomes on mode b out , which corresponds to tracing over this mode. This leads to 


+oo 


+oo 


i$>out = i$> fe , i$>* = £\ 

k =0 n=k 

The photon statistics are then given by the probability distribution 


t n k \toL\ n — k) . 


( 10 ) 


+oo 


ip(m) = y^r 2k k\ 


k=0 


+oo 

E< 

n=k 


t n k (m\ta; n — k) 


In order to express ip(m) more explicitly, we recur to the relation 

D(a)(o t ) n |0) = (a t —a*)”£>(a)|0), 


( 11 ) 


( 12 ) 


which is verified by virtue of the unitarity of the displacement operator D, along with the relation W (a)a) D (a) = a) + a*. 
Consequently, the inner product (m\a; n) which appears in Eq.(l 1) can be written as 


{m\a; n) = (m| £ (") (o t )"(-a*)'*-'‘|a) = 

re-0 V'V 


= E 


«=o 

n \ \fm\ 


(13) 




(^_^n-K\ a yn+n-2K e i(m-n)4) e ~\a \ 2 /2 ^ 


where a = \a\e 1 ^. 

Finally, we have to take into account the response of the time-multiplexed detector (TMD). This is characterised by its non¬ 
unit detection efficiency 77 , which we model as a beam splitter with reflectivity 77 , and draws its photon-number resolution 
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from temporal binning. These two aspects - loss and binned photon counting - are accounted for by two matrices: L(rj) 
models photon loss, while C describes the binning operation [10, 21, 46]. If we denote by f = {<^(0), <p(l), <p(2 ),...} the 
vector of photon statistics and by p = {p(0),p(l),p(2 ),...} the vector of counting statistics, the following relation holds: 

P = C- L(rj) • f. (14) 

This single-mode expression can be easily generalised to our two-mode layout. The joint photon statistics for two WFH 
detectors with settings {r, £, a} and {r', a observing the state c n,n' \ n )\ n ')i are given by 


+oo +oo 


4>(m,m / ) = k\k'\r 2k (r') 2k 


k =0 fc'=0 


+oo +oo 

EE' 

n=k n'=k' 


t n k (t') n k (m\ta; n — k)(m'\t'a'; n' — k') 


(15) 


Loss on each mode - described by the efficiencies rft and r] 2 , respectively - and temporal binning can be taken into account 
by casting the probabilities T>(m, m!) into a matrix F, which is related to the matrix of the joint detection probabilities P 
through [46] 


E = Ci-L(m)-E-L(rti) T -Cj. 


(16) 


Supplementary Note 2: Theoretical model of the experiment 

We adopt a collinear geometry for our WFH detectors in order to simplify the experimental scheme and ensure better 
passive phase stability (see Supplementary Fig. 1 and Supplementary Fig. 3). The relative phase between signal and local 
oscillator is varied by means of birefringent elements arranged in a geometric phase rotator (GPR). However, imperfections 
in the individual wave plates and calibration inaccuracies can determine an undesired modulation of the reflectivity associated 
to the polarising beam splitter in each WFH detector. We thus perform a calibration of both GPRs in the presence of the local 
oscillator only. A typical plot of the dependence of counts on the setting of the GPR is shown in Supplementary Fig. 2. We 
account for this spurious effect by expressing the reflectivity as r(0) = roy/l + v cos(40 + 0 O )> where ro = 0.5 is the ideal 
reflectivity, 0 is the setting of the half-wave plate in the GPR, 6o is an offset angle, and the modulation depth v is expected to 
be of the order of 0.1 for small deviations from the perfect case. This relation for r(0) is included in Eq.(16) for the expected 
joint detection probabilities. 

Modelling a split single-photon state (SSPS). We want to compare our experimentally generated split single photon to a 
theoretical signal. It is reasonable to assume that the SSPS exhibits some degree of mixedness [26], therefore we take into 
account imperfect single-photon generation by including contributions from the vacuum and higher-order terms from the 
downconversion process (given the brightness of our source) in the expression of our state: 

Po = w 0 |0>(0| + wi|l)(l| + (1 ~ w 0 - wi)|2)(2|. (17) 

This state produces an entangled resource when it is split on a symmetric beam splitter. The corresponding counting statistics 
read 


E = w 0 Po + w 1 Pi + (1 - Wo - wi)P2, (18) 

where the matrix P describes the statistics associated with the entangled state obtained from Fock layer | k ). 

The coefficients wo and wi, along with the detection efficiencies t]a and tjb , are the key parameters which appear in 
our model. The latter are determined from the counting statistics of the experimental input state when no local oscillator 
impinges on the PBS, following the Klyshko calibration method. In order to compute wq and wi, we suitably invert Eq.(16) 
and obtain the photon statistics matrix of the source, which we label F PDC ; we thus estimate wo = F PDC (0, 0) and w\ = 
Z PDC (1,0)+Z PDC (0,1). 

In Supplementary Fig. 4 we show the theoretical curves for the joint counting statistics matrix P up to two detection events 
on each output mode. These plots inform us about the main features of the physical phenomenon which we probe with our 
WFH detectors: coherent oscillations across Fock layers are present in all terms but those where at least one side detects 
no photons. Moreover, the oscillations are shown to be in phase with one another. As a further investigation of the role of 
loss in our scheme (which we have already taken into account by using the experimental detection efficiencies in our model), 
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Supplementary Fig. 4 includes three distinct detection scenarios given by different values of 77 a and 77 b ; poor efficiency of 
detection appears to mainly affect the visibility of the studied nonclassical oscillations. 

Modelling a two-mode squeezed state (TMSS). The expression for an ideal two-mode squeezed vacuum state reads 

+00 

|^)tMSS = \/l l^l 2 E ^ n \ n )^\ n )Bi (19) 

n=0 

where |A| is the real squeezing parameter determined by the brightness of the parametric downconversion source. However, 
a more realistic expression for the experimentally produced state is 

p = (1 - p)|^)tmss(^| +p|0,1) (0,1|, (20) 

where the term |0,1)(0,1| is a first-order approximation of noise in a squeezed thermal state. Indeed, the presence of such 
contribution is suggested by the experimental g^ for the source, whose values highlight the two following features for our 
squeezer: an asymmetry across the two modes, and the presence of unwanted frequency modes in the generated two-mode 
squeezed state. Both effects are captured by the additional term appearing in Eq.(20). 

In this case, the relevant quantities for the theoretical model are given by the squeezing parameter and the weight p. The 
former can be computed from the photon statistics of the downconversion source; manual inspection of the experimental and 
theoretical plots suggests the value p = 0.04 for the coefficient in the expression of p. 

The theoretical curves for the joint counting statistics matrix P (up to two detection events) for the TMSS are presented in 
Supplementary Fig. 5. Similarly to the case of a split single photon, the model predicts coherent, in-phase oscillations across 
the Fock layers for all terms but those where at least one side detects vacuum. Once again, the effect of loss is a reduction of 
the visibility of such oscillations. 

As already pointed out in the main text, we remark that our theoretical models the SSPS and the TMSS share some common 
limitations. Notably, imperfect mode matching affects the interference between the probed state and the classical reference. 
For this reason, this aspect is likely to act adversely on the joint detection statistics too. Nevertheless, our descriptions have 
the advantage of relying on very few parameters, most of which are directly estimated from the collected data. Each model 
agrees with the data for all considered terms in the joint counting statistics simultaneously, which in turn means that we can 
harness the behaviour of all multi-photon components in a given output state with one single description and set of parameters. 


Supplementary Note 3: Nonlocality of a two-mode squeezed state detected with weak-field homodyne 

In this section we provide details on the application of the Collins-Gisin-Linden-Massar-Popescu (CGLMP) inequali¬ 
ties [33] to WFH. Following the notation introduced in Supplementary Fig. 6, we illustrate the case in which the number of 
detected photons on each side is M = 2, i.e., we restrict our study to local subsystems with dimension D = M + 1 = 3. The 
generalisation to higher photon numbers proceeds in a similar fashion. Fet us consider the two-mode quantum state |4 /)tmss> 
distributed across modes a|, a\, and the two local oscillators | a) and \f3) on modes b[ and b\ respectively. As we focus on 
the detection of low photon numbers (M = 2), we write the input state as 

Win = (b\) 2 (bl) 2 + y(al) 2 (4) 2 + Ao^a^^j |0) s |0)i|0)i|0) 2 , (21) 

where a = \a\e l ^ a , /3 = \/3\e l( ^ b and A = |A|e^. 

The BS of each WFH detector operates the following transformation from the input modes into output modes {d\ , e\ , d\ ^2 } : 

a| = t\d\ + ri, 
b\ = nd\ - tie\ , 

4 = + r|4 , 

„ 4 = r 2 4 - *2 4 > 


< 


( 22 ) 
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where we assume that tj E M and ry = —i\rj \ for j = 1,2. 

These relations lead to the expression for the output state |\k) out in the Fock basis (up to a normalisation factor): 

Wout = (—y-r?r 2 + y*i*2 + Aa/?tit 2 rir 2 ) |2200) + + yW) 2 ^) 2 + Xaptit 2 rlr^j |0022) + 

+ + y^(r 2 ) 2 - XaPtit 2 rir^j |2002) + (^y-^ r 2 + y(^) 2 t 2 - Xaptit 2 r\r^j |0220) 

+ [a: 2 /? 2 iif 2 rir 2 + 2X 2 t 1 t 2 r* l rl - Aa/?(i 2 |ri| 2 - |ri| 2 |r 2 | 2 + i 2 |r 2 | 2 - t\t 2 2 )] |1111) + 

2 o2 

+ -—^-hr\r 2 + X 2 t\t 2 r 2 - Aa/S^if^i - firi|r 2 | 2 ) |2101) + 

2 o2 

+ y-iifir 2 + A 2 fif 2 ri - XaP(t\t 2 r 2 - t 2 \r\\ 2 r 2 ) 11210) + 

2 /q2 

+ tit 2 r 2 + A 2 t 2 (r^) 2 r2 - Aa/3(tiri|r 2 | 2 - titirj) |0121)+ 

r ry 2 ft 2 1 

+-tit 2 ri + A 2 tiri(r2) 2 + \aP(t\t 2 rl - t 2 \ri\ 2 rl) |1012). 


(23) 


Here we use the convention |0000) = |0)i t |0) 2t |0)i r |0) 2r , and similarly for the other number states (‘t’ and ‘r’ denoting 
transmitted and reflected output modes, respectively). 

In the work of Collins and collaborators, the local realistic bound is calculated on a linear combination of detection prob¬ 
abilities associated to local measurement settings - in our case these are the complex parameters a and /3. As introduced 
in [33], we adopt the shorthand notation P(Ai = Bj + e) to indicate the probability of an event for which the detection 
outcomes Ai and Bj differ by e = — 1, 0,1 (where { — 1, 0,1} and {0,1, 2} are congruent modulo 3). Therefore, the CGLMP 
inequality adapted to our layout reads 


\h\< 2, 


where 


h =[P(A i = Bi) + P(B 1 =A 2 + 1)+ P{A 2 = B 2 ) + P(B 2 = Ai)]- 

- [P(i4i = B 1 - 1) + P(Bi = A 2 ) + P(A 2 =B 2 - 1) + P(P 2 = Ai - 1)]. 


According to the notation introduced above, each term in Eq.(24) is computed as follows: 

P(Ai = B\) = (Poo )( a i! fix') + (Pn )(o:i> Pi) + (P22 )(o ; ij Pi), 
P(Bi = A 2 + 1) = ( P10 )(o: 2 , Pi) + ( P21 }(o: 2 , Pi) + (Po 2 )(a 2 , Pi), 
P(A 2 = B 2 ) = (P 0 o)(a 2 , P2) + (At )( a 2, P2) + (P22 ){ot 2 , p 2 ), 
P(B 2 = A\) = (P„o>(ai,&) + (At )(oti,p 2 ) + ( P22 ){oii,p 2 ), 
P(Ai = Bi — 1 ) = (P10 )( a i,Pi) + (P21 )ipi,Pi) + (Pq 2 )(ai,Pi), 
P(Bi = A 2 ) = ( Poo )(<a 2 , Pi) + (Pn )(<a 2 , Pi) + (P22 )(°-2, Pi), 
P(A 2 = B 2 — 1 ) = (P10 )(ot 2 ,p 2 ) + (P 2 1 )(ai 2 ,p 2 ) + (P02 )(u 2 ,p 2 ), 
P(B 2 = Ai - 1) = (P01 )(«i)# 2 ) + (P12 )(cki ,/?2) + (P 2 o)(ai,p 2 ), 
where we conveniently relabelled the four-photon output states. Namely, 


12200) <—► |00) ->■ Poo , 
|0022) <—► |22) P 22 , 

|0121) <—► |21) ->■ P 21 , 
|2101) <—► |01) -)> Pot , 
|0220) <—► 120) -»• P 20 . 


11111) <—> 111) —>■ P 11 , 
|1210) <—► 110) -> P 10 , 

12002) <—► 102) -> P 02 , 

|1012) <—► 112) -> P 12 , 


(24) 


(25) 
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This operation amounts to retaining the number of detected photons on the two reflected output modes as the two significant 
digits for the labels. The expectation values (•) of the projection operators {Pkk'}> fc, fc' = 0,1, 2, appear in Eq.(25). These 
are calculated on the output state |\P) out , i.e., (Pkk' ) = {^\Pkk' |4/) ou t f° r all k,k' considered. The quantities {a/,/?//}, 
1,1' = 1,2, refer to physical measurement settings, which in our case are identified by the complex amplitude of each local 
oscillator. 

In Supplementary Fig. 7 we show the oscillations of the terms {(Pkk' )}, k, k' = 0,1, 2, which appear in Eq.(24) through 
the detection probabilities P(Ai = Bj + e), as the phase of one local oscillator is varied. Once again, we note that this 
study assumes an ideal setup, from the generation of pure two-mode squeezed states to perfect WFH detectors. The most 
remarkable feature in Supplementary Fig. 7 is the presence of two different oscillation periods for {(Pkk ')}, the period being 
determined by the specific Fock layer observed. 

The presence of different oscillation periods associated to distinct photon numbers has been observed when a single-mode 
squeezed state is interfered with a coherent state of comparable energy [47]. However, the presence of such super-resolved 
fringes does not imply a quantum enhancement in phase estimation. In fact, while super-resolution can be obtained by purely 
classical means [48, 49], super-sensitivity can only be achieved with quantum resources. An analysis in this direction needs to 
account for all resources and all post-selection probabilities. The usefulness of the scheme discussed in this work for quantum 
metrology is currently under investigation. 

Finally, we comment on the role of detection efficiency in the framework of the generalised Bell inequalities we have just 
presented. As shown in Supplementary Fig. 4 and Supplementary Fig. 5, lower detection efficiencies decrease the visibility of 
the oscillations and artificially enhance the low-photon components. The effect is clearly more detrimental than the induction 
of a detection loophole, as in standard Bell tests performed with entangled photons. 

Further, this effect is more pronounced for higher Fock layers. In fact, the first-order contribution to noise in Fock layer M 
comes from the higher layer M + 1. Each term in Eq.(24) may then be written as 

P oc r/ M P M + (M+1)t? m (1 - rf)v M +i (26) 

where Pm is the correlation in the ideal case r]=l , and vm+i is a noise term coming from Fock layer M + 1. The latter is 
due to the loss of a single photon, which one could otherwise reject by post-selection. We thus see that the relative weight 
becomes more important as M increases. Therefore, the achievable violation not only decreases with M, but its resilience to 
loss also becomes weaker. 
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